# include "mempool.hpp"
# include "taskpool.hpp"
# include <iostream>

# include <new>
# include <vector>
# include <pthread.h>

# include <ctime>

# include "linalg.hpp"

# include "geometry3d.hpp"
# include "geometry2d.hpp"
# include "stdio.hpp"

# include <algorithm>
# include <iterator>

# include <cstdlib>


namespace test {


  bool testMatrix() {
    using linalg::Matrix;
    using linalg::Array;

    Matrix m(2, 2);
    m(0, 0) = 1;
    m(0, 1) = 2;
    m(1, 0) = 3;
    m(1, 1) = 4;

    double const* const data = m.data();
    std::copy(data, data+4, std::ostream_iterator<double>(std::cout, " "));
    std::cout << "\n";

    Array r1 = m.row(0), r2 = m.row(1);

    std::cout << r1(0) << " " << r1(1) << std::endl;
    std::cout << r2(0) << " " << r2(1) << std::endl;


    Array c1 = m.column(0), c2 = m.column(1);

    std::cout << c1(0) << " " << c1(1) << std::endl;
    std::cout << c2(0) << " " << c2(1) << std::endl;

    return true;
  }

  bool testSVD() {
    using linalg::Matrix;
    using linalg::svd;

    Matrix m(2, 2);
    m(0, 0) = 1;
    m(0, 1) = 2;
    m(1, 0) = 3;
    m(1, 1) = 4;

    Matrix u, s, vt;
    svd(m, u, s, vt);

    return true;
  }

  bool test3D() {
    using geometry::d3::Point;
    using geometry::d3::Vector;

    Vector v = {1, 2, 3};
    Vector u = {-1, -2, -3};

    Vector t = u+v;
    std::cout << t.x() << std::endl;

    return true;
  }

  bool testFitting() {
    using geometry::d3::Point;
    using geometry::d3::CoordinateSystem;
    using geometry::d3::computeFittingPlane;

    std::vector<Point> ps;
    for (size_t i = 0; i < 40; ++i)
      ps.push_back(Point(drand48(), drand48(), drand48()));

    auto p = computeFittingPlane(ps);
    std::cout << p.origin << std::endl;

    return true;
  }

  void printCurve(std::vector<geometry::d2::Point> const& ps) {
    for (geometry::d2::Point const& p : ps)
      std::cout << p.x() << " " << p.y() << std::endl;
  }

  bool testCurveSimplification() {
    using geometry::d2::Point;
    using geometry::d2::simplifyCurve;

    std::vector<Point> c;
    for (int i = 0; i < 10; ++i)
      c.push_back({(double)i, i+(drand48() - 0.5)});

    printCurve(c);
    std::cout << "\n\n\n";

    std::vector<Point> const sc = simplifyCurve(c, 1.e-6);
    printCurve(sc);

    return true;
  }

  bool testSolve() {
    using linalg::Matrix;
    using linalg::solve;

    Matrix A(2, 2);
    A(0, 0) = 1;
    A(0, 1) = 5;
    A(1, 0) = 7;
    A(1, 1) = 17;

    Matrix b(2, 2);
    b(0, 0) = 0;
    b(1, 0) = 1;

    b(0, 1) = 1;
    b(1, 1) = 0;

    A.LUD();
    solve(A, b);

    std::cout << b(0, 0) << " " << b(1, 0) << std::endl;
    std::cout << b(0, 1) << " " << b(1, 1) << std::endl;

    return true;
  }

}


# include <type_traits>

int main() {

  test::testMatrix();

//    std::cout << "point 3d " << sizeof(geometry::d3::Point) << std::endl;
//    std::cout << "point 2d " << sizeof(geometry::d2::Point) << std::endl;

//    std::cout << "vector 3d " << sizeof(geometry::d3::Vector) << std::endl;
//    std::cout << "vector 2d " << sizeof(geometry::d2::Vector) << std::endl;

//    std::cout << "double x 3 " << sizeof(double) * 3 << std::endl;
//    std::cout << "double x 2"  << sizeof(double) * 2 << std::endl;

//    test::testCurveSimplification();

//    test::testSolve();

    return 0;

}


